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SUMMARY 


A comparative study of two different incompressible Navier-Stokes algorithms for 
solving an unsteady, incompressible, internal flow problem is performed. The first algo- 
rithm uses an artificial compressibility method coupled with upwind differencing and a line 
relaxation scheme. The second algorithm uses a fractional step method with a staggered 
grid, finite volume approach. Unsteady, viscous, incompressible, internal flow through a 
channel with a constriction is computed using the first algorithm. A grid resolution study 
and parameter studies on the artificial compressibility coefficient and the maximum allow- 
able residual of the continuity equation are performed. The periodicity of the solution is 
examined and several periodic data sets are generated using the first algorithm. These 
computational results are compared with previously published results computed using the 
second algorithm and experimental data. 


INTRODUCTION 


In the field of computational fluid dynamics, it is challenging to obtain accurate re- 
sults in a computationally efficient manner. A wide variety of approaches can be employed 
when solving a specific problem. This is particularly true of unsteady, incompressible meth- 
ods; there are no time derivatives in the continuity equation, yet the velocity field must be 
divergence free at every time step. Manipulation of the governing equations must occur in 
order to formulate the problem. The current study compares two different incompressible 
Navier-Stokes algorithms to solve an unsteady, incompressible, internal flow problem. The 
first approach uses the artificial compressibility method introduced by Chorin (ref. 1) cou- 
pled with upwind differencing and a line relaxation scheme. The second approach uses a 
fractional step method (refs. 2 and 3) with a staggered grid, finite volume method. These 
two approaches have been widely used in recent work. As an example, in the 1991 10th 
AIAA Computational Fluid Dynamics Conference, a number of incompressible flow papers 
which used these methods were presented. A number of authors (refs. 4-9) used the ai Un- 
cial compressibility method in their work presented at this conference. An equal number 
of authors (refs. 8-13) used the fractional step method as a basis for their formulations. 



Both methods have been used to compute the unsteady, viscous, incompressible, 
internal flow through a channel with a constriction. This specific problem was chosen 
for two reasons: one is the availability of experimental data for the problem; another is 
that this problem is geometrically simple but produces complex structures in the flow 
field for analysis. The physical applications of this problem include the area of bio-fluids 
research. Pulsatile flow through a channel with an obstruction models the phenomena 
of arteriosclerosis closely. This study provides not only the opportunity to perform code 
validation and comparison with experimental data, but generates results which may aid in 
a better understanding of the fluid dynamics associated with arteriosclerosis. 

The current work includes a grid resolution study to determine an appropriate grid 
for the flow field analysis. A parameter study was performed on the artificial compress- 
ibility coefficient and its effect on solution convergence and accuracy. The solution was 
examined to determine the number of physical time steps needed to obtain a periodic 
solution. Solution accuracy and the effect of varying the maximum allowable tolerance 
for the divergence of velocity was investigated. Finally, the work was compared with the 
computational data of Rosenfeld (ref. 14) and the experimental results of Park (ref. 15). 

The authors wish to thank Dr. Moshe Rosenfeld of MCAT Institute, San Jose, for 
his generous contributions of time and experimental results to this paper. 


METHODS 


The incompressible, unsteady, Navier-Stokes equations are known to be particularly 
difficult to solve. One of the reasons for this difficulty is the absence of a time derivative 
of pressure in the governing equations. Since the speed of sound in a truly incompressible 
fluid is infinite, the entire pressure field is affected instantaneously by a disturbance at 
any one point in the flow domain. This requires that the numerical algorithm propagate 
information through the entire flow domain during one discrete time step. Global coupling 
of the discretized system is implied, thus some type of iterative scheme is usually required 
to solve the equations in time. The elliptical nature of these equations causes the generation 
of time-accurate solutions to be computationally expensive. 

Two different computational methods were employed to model the fluid dynamics of 
unsteady, incompressible, pulsatile channel flow: (1) an artificial compressibility method, 
and (2) a fractional step method. The fractional step data used for comparison in this 
study was generated separately by Rosenfeld (ref. 14), by the fractional step code, INS3D- 
FS (ref. 16). The INS3D-FS code utilizes an integral formulation of the Navier-Stokes 
equations and finite volume discretization to solve both steady-state and time- dependent 
flows. In this algorithm, the momentum equation is marched forward using the last known 
pressures to solve for the velocity field. The new velocity field at this point does not 
necessarily satisfy the continuity equation. A Poisson equation in pressure is then solved 
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to obtain the corrections to the velocity field such that the continuity equation will be 
satisfied; thus, the solution is advanced one time step. In practice, the Poisson equation is 
iterated until the maximum residual is less than 10 8 , where the residual is the divergence 
of velocity multiplied by the volume of the computational cell. A staggered mesh approach 
is used with the INS3D-FS code, which removes the need for artificial dissipation. 

The code used by the authors to generate solutions was the INS2D code developed by 
Rogers (refs 17 and 18). This code solves the incompressible Navier-Stokes equations in 
two dimensional generalized coordinates for both steady-state and time varying flows. The 
equations are formulated into a mixed set of hyperbolic and parabolic partial differentia 
equations using the method of artificial compressibility. In this method, a pseudo-time 
derivative of pressure is added to the continuity equation. 


dp! dr — —flVlf 


where p = pressure, r = pseudo-time, if = velocity, and /? is the artificial compressibility 
constant which controls the propagation of the artificial pressure waves. The convec- 
tive terms are differenced using an upwind biased flux-difference splitting method. The 
discretized equations are solved using a line-relaxation scheme. This method uses sub- 
iterations in pseudo-time on the global system of equations as the mechanism to drive the 
divergence of velocity toward zero. In practice, the sub-iterations are terminated once the 
maximum divergence of velocity in the continuity equation is reduced below a specified 
constant, £ cont . This constant can be much larger than machine zero; for instance, accurate 
solutions have been obtained (refs. 17 and 18) using values of e C ont — 10 


PROBLEM DEFINITION 


The physical problem examined in this study was pulsatile flow through a channel 
with a constriction. The computational geometry is consistent with the experimental set- 
up of Park (ref. 15), and is shown in figure 1. The channel height is h and has been 
normalized to unity. The height of the constriction is given by a = 0.57. This is the 
distance from the top wall of the channel to the lowest point in the constriction. The 
length of the channel upstream of the constriction is given by L u = 7. The length of 
the channel in which the constriction occurs is given by L c = 4.66, and the downstream 
portion of the channel is given by L d = 15.34. The total length of the channel is given by 
Ltotal and is equal to L u + L c + Ld = 27.0. 

The computational boundary conditions were implemented as follows. For the upper 
and lower walls of the channel, no slip conditions were specified. The inflow boundary 
to the experimental set-up occurred at 100 channel heights upstream of the constriction. 
Since it would be computationally expensive to compute the flow for this entire domain, 
the computational inflow boundary was placed at seven channel heights upstream of the 
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Figure 1. Physical description of computational geometry. 


constriction. Previous work was done by Rosenfeld (ref. 14) to determine the appropriate 
inflow boundary conditions to match the experimental setup, and this information was 
used in the current work. A parabolic profile was imposed at the inflow boundary and the 
average inflow velocity was set to match the mass flow from the experimental setup. The 
average inflow velocity is periodic and given by the shape function in figure 2. 



Figure 2. Average inflow velocity magnitude versus time for one period. 
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This shape is defined analytically by 

17(f) = U S 0 < t/T < 1/2 

17(f) = U s -U p sin (27rf/T) 1/2 < t/T < 1 

where U s is the non-dimensionalized steady component of the average inflow velocity, I p 
is the pulsatile component and T represents the period. This waveform was chosen for the 
inflow velocity because it approximates the pulses in velocity that blood encounters during 
diastole and systole in mammalian circulatory systems and provides realistic conditions 
for the modeling of arteriosclerosis. 


Associated with this problem, there are several relevant relationships between the 
physical parameters which govern this flow. 


v 

Re s 

~T~ 

h 

U^f 

= hfp_ 
p StRe.. 


fp ~ 


St = 


Re s is the Reynolds number for the flow, U s is the steady component of the average inflow 
velocity, h is the non-dimensionalized height of the channel, and u is the kinematic viscosity 
of the fluid. The frequency of the period T is given by f p . The Strouhal number is given by 
St, where h is the channel height, and U p is the pulsatile component of the average inflow 
velocity. These physical parameters were specified to correspond with the computational 
work done by Rosenfeld (ref. 14), and the experimental work done by Park (ref. 15). 


The outflow boundary conditions specified in the work by Rosenfeld (ref. 14) were of 
the non-reflecting type with the streamwise derivative of velocity set to zero, du/dx = 0. 
In the present work, a non-reflecting boundary condition was used based on the method 
of characteristics to solve for the velocity; the static pressure was specified to be constant. 


To compare the computed results from the two algorithms on an equal basis, the 
grids used to model this geometry matched the grids used by Rosenfeld (ref. 14). A coarse 
mesh was generated with 97 X 21 grid points and clustering near the upper and lower walls 
and in the regions where large gradients were expected, resulting in 70 percent of the grid 
points being clustered in the vortical region downstream of the constriction. From this 
mesh, medium, fine, and ultra-fine grids were generated by doubling the points in each of 
the i and j directions, as defined in figure 1, for each level of refinement. This resulted in 
four grids with an increasing number of grid points: a coarse grid with 97 x 21 grid points, 
a medium grid with 193 x 41 grid points, a fine grid with 289 x 61 grid points, and an 
ultra-fine grid with 385 x 81 grid points. Figure 3 shows part of the coarse mesh near the 
constriction. 
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A baseline case was first calculated from which further parametric studies were per- 
formed. This case shall be referred to as Case I, and the constants specified are: the 
Reynolds number Re s = 131.9, the depth of the constriction a/h = 0.57, the Strouhal 
number St = 0.42, the pulsatile component of the average inflow velocity U p — 1.217, 
and the period T = 1.957 seconds in physical time. 



Figure 3. 97 x 21 coarse mesh in the area of the constriction. 


COMPUTED RESULTS 


First, a grid resolution study was performed using all four grids and the parameters 
of Case I. To obtain a starting solution for each of the grids, a converged steady-state 
solution was obtained. The calculation was restarted with time-varying pulsatile inflow 
from the steady-state solution and run out for four periods with a value of (3 = 50, 000 
arid At = 0.03914. Grid convergence was analyzed using two physical parameters: the 
distribution of skin friction on the bottom wall of the channel, and the static pressure at 
the centerline of the channel. As these are laminar calculations, the skin friction coefficient, 
C/, is normalized by the Reynolds number and is given by 

n — Tw 

’ 1/2 P U? 

where r w is the shear stress at the wall and p is the density. The skin friction coefficient 
and pressure at the centerline of the channel are shown at the beginning of the period, 
t/T = 0.0, for all four cases shown in figures 4(a)-(b) and 5(a)-(b). From the figure, it is 
apparent that the solutions are close to being grid independent for the two finest meshes, 
as the solution from the 289 x 61 mesh is close to the solution from the 385 x 81 mesh. 
However, 30 percent less CPU time is needed to generate a solution for a 289 x 61 mesh 
than for the 385 x 81 mesh. This geometry produces a solution with many strong flow 
gradients. Flow is produced which is continuously reversing throughout the channel which 
creates strong streamwise pressure gradients and strong vorticity gradients near the walls. 
A grid with relatively fine spacing is needed to adequately capture all of the important 
flow physics which are inherent in this problem. As a result, the 289 x 61 grid was chosen 
to use in the remainder of the calculations in the study. Rosenfeld (ref. 14) performed 
a grid resolution study and determined that a nearly grid independent solution was also 
obtained with the fractional step method for the 289 x 61 grid. 

An in-depth examination of the skin-friction and pressure data, shown in figures 4(b) 
and 5(b), reveals several “spikes” in the INS2D solution near the beginning and end of 
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Figure 4(a). Effect of grid resolution on skin friction at the lower wall. 



the constriction. In these figures, the magnitude of spikes appears to diminish as the grids 
become more refined. The magnitude of the spikes in the solutions was also affected by 
the choice of ft, the artificial compressibility coefficient, used in running the INS2D code. 
This dependence on ft is discussed in the following section. 

A parametric study on the effects of different values of ft on the speed of solution 
convergence was performed. The values of ft used ranged from 50,000 to 0.1 by an order of 
magnitude for each case. When running the INS 2D code, the use of different values of ft 
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Figure 5(a). Effect of grid resolution on centerline pressure. 



Figure 5(b). Close-up of figure 5(a), in area of constriction. 


will generate basically the same solution with minor differences, but will greatly affect the 
speed of convergence and therefore the amount of CPU time needed to obtain a solution. 
This occurs because the number of sub-iterations that are calculated in pseudo-time in 
the unsteady mode axe either increased or decreased depending on how a specific value of 
[j propagates information within a specific geometric configuration. Figure 6 shows the 
effect of different values of (3 on CPU time; the data from these solutions is shown in fig- 
ures 7(a)-(b) and 8(a)-(b). The use of a large (3 greatly enhances the speed at which 
the unsteady solution converges for this geometry, as shown in figure 6. However, in 
figures 7(a)-(b) and 8(a)-(b), which show the skin-friction coefficient at the lower wall and 
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Figure 6. Effect of 0 on computing time needed for first half of period. 
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Figure 7(a). Effect of 0 on skin friction at the lower wall at t/T = 0.5. 

pressure at the centerline of the channel at t/T = 0.5, the relative smoothness of the 
solution is adversely affected by the use of the large 0. The pressure appears to be more 
adversely affected by the use of a large 0 than the skin-friction. Since the pressure term 




in the continuity equation is scaled by ft, small spikes in the velocity field produce large 
spikes in the pressure field. A value of ft = 10 3 was chosen for use in the remainder of 
the study. This value produced reasonably fast convergence to a periodic solution and 
minimized the magnitude of the spikes in the solution. 



Figure 7(b). Close-up of figure 7(a), in the area of the constriction. 


Several alternatives can be utilized to reduce the magnitude of the spikes in the 
solutions of the large ft cases. Since the grid resolution study, shown in figures 4(b) 
and 5(b), shows that the use of a finer grid reduces the magnitude of the spikes, refining 
the grid in the area of the constriction could produce a smoother solution with the use 
of a larger ft. This could be done without increasing the number of points in the grid by 
clustering the existing points in the constriction region. Also, smoothing of the grid was 
performed by fitting a cubic spline to the grid lines lengthwise at the beginning and end 
of the constriction. This smoothing did reduce the magnitude of the spikes for a case of 
ft = 50, 000, as shown in figures 9(a)-(b). However, since this is a validation study with the 
results of Rosenfeld (ref. 14), the remainder of this study was performed with the exact 
grids used in the work done by Rosenfeld (ref. 14) and not the smoothed grids. 

The use a TVD limiting scheme is another method which may be employed to re- 
duce the magnitude of the spikes when running the INS2D code with a large ft. Using a 
TVD limiting scheme with the upwind differencing may add enough dissipation to smooth 
out these spikes, yet may also reduce the accuracy of the solution. Thus, the preferred 
correction for this case would involve refining the affected region of the grid. 

The issues of periodicity and convergence of the solution were also examined in this 
study. When running a time dependent case with the INS2D code, a convergence critera 
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Figure 8(a). Effect of 0 on centerline pressure at t/T = 0.5. 



Figure 8(b). Close-up of figure 8(a), in the area of the constriction. 

must be specified for the sub-iterations. The parameter e conU the residual of the continuity 
equation, is used as the tolerance for the maximum divergence of velocity in the entire 
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Figure 9(a). Smoothed versus non-smoothed grid for (3 = 50,000. 



Figure 9(b). Close-up of figure 9(a), in the area of the constriction. 

flow field; once the tolerance is reached, the sub-iterations are terminated and the code 
continues on to the next time level. The maximum divergence of velocity at each grid point 
was determined using a third order differencing scheme. A parameter study was performed 
to determine the appropriate value for e CO nt which would generate an accurate periodic 
solution for the least amount of computer time. The value of e CO nt was varied by an order 
of magnitude for four cases from 10 -1 to 10 -4 . Solutions were generated for each value of 
£cont for 10 consecutive periods. The solutions were checked for periodicity by taking the 
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root mean square (RMS) of the difference over each grid point of the dependent variables 
between subsequent periods. A direct correlation exists between the decrease m these R 
values and the solution periodicity. Figure 10 shows the effect of e cont on the penodici y 
of the solutions. As shown in figure 10, running the INS2D code with a value of £ cont 
of 10~ 3 produces the solution with the smallest RMS value consistently throughout the 
ten periods. Figure 11 shows the same data plotted against CPU time instead of periods. 
Each order of magnitude decrease in e cont produced a factor-of-two increase m the number 
of sub-iterations within one time step, and therefore increased the amount of CPU time 
needed to obtain a solution. Values of £ cont of 10' 3 and 10~ 4 generated the more periodic 
solutions, but a value of e cont = 1(T 3 produces a smaller RMS error than e cont i - 10 > 

this translates to faster convergence to periodicity, while using approximately 40 percen 
less CPU time. The value of £ c0 „< of 10 -3 was used in the remainder of the study. 



Figure 10. Effect of e con t on periodicity: root mean square of change in dependent variables 
versus period. 


Once the basic operational parameters had been established for running the INS2D 
code with this configuration, two cases were run to generate data which were used to 
compare with Rosenfeld’s (ref. 14) computational data and Park’s (ref. 15) experimental 
data. The geometry remained the same for the two cases, but the inflow conditions vary 
and axe given in table 1. 
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Figure 11. Effect of e con t on periodicity: root mean square of change in dependent variables 
versus CPU time. 


Table 1. Parameters used for comparison study 


Parameter 

Case I 

Case II 

Rz s 

131.9 

138.3 

ajh 

0.57 

0.57 

T 

1.957 

1.368 

St 

0.42 

0.42 

u 3 

1.5 

1.5 

U p 

1.217 

1.740 


A converged periodic solution was generated for both of these cases. Data was saved 
at each time step in the period for comparison with Rosenfeld (ref. 14) and Park (ref. 15). 
Figures 12(a)-(d) show the coefficient of skin friction on the lower wall for both the current 
results using the INS2D code and Rosenfeld’s (ref. 14) work using the INS3D-FS code, at 
various times during the period for Case I. A comparison of the computational data shows 
that the solutions appear to match with respect to location of skin friction maxima and 
minima. However, the data generated using the INS3D-FS code has a larger fluctuation 
in skin friction values in the region just behind the constriction than the data generated 
using the INS2D code. This correlates to stronger vorticity in these regions in the solution 
generated by the INS3D-FS code. 


14 





Figure 12(a). Comparison of INS2D and INS3D-FS results: computed skin friction at the 
lower wall at t/T = 0.24. 



Figure 12(b). Comparison of INS2D and INS3D-FS results: computed skin friction at the 
lower wall at t/T - 0.5. 


The INS2D code was also used to compute instantaneous streamlines for Case I. 
These streamlines are compared with streamlines generated using the INS3D-FS code for 
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0.0008 



Figure 12(c). Comparison of INS2D and INS3D-FS results: computed skin friction at the 
lower wall at t/T — 0.74. 



Figure 12(d). Comparison of INS2D and INS3D-FS results: computed skin friction at the 
lower wall at t/T = 1.0. 

Case I and are shown at 0.1T intervals throughout the period in figures 13 and 14. The 
streamline plots illustrate the same phenomena as the skin friction plots, that the results 
generated by the INS3D-FS code have larger vortical structures behind the constriction 
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Figure 13. INS2D results: streamline data at 0.1 T intervals from 0.1T to 1.0T. 

throughout the entire period, even during the pulse in the inflow average velocity, as shown 
in frame 7 at 0.7T in figures 13 and 14, when most of the vortical structures are washed 

out. 

Another easily quantified physical result was the location of the center of the 
B-vortex The B-vortex is defined to be the vortex along the bottom wall of the channel; 
it grows immediately behind the end of the constriction and is shed downstream. The lo- 
cation of this vortex was examined as function of time and the results were compared wit] 
the INS3D-FS code results and the experimental results for both Case I and Case II. 1 he 
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Figure 14. INS3D-FS results: streamline data at 0.1T intervals from 0.1T to 1.0T. 


location of the center of the B-vortex can be found by plotting the instantaneous stream- 
lines with a fine increment of the contour lines. This is analogous to the measurement of 
the experimental vortex location, which was measured using flow visualization pictures. 
Figures 15 and 16 show the results given by Rosenfeld (ref. 14), the computed results from 
this study, and the experimental results. There is good agreement between the computed 
results by the INS3D-FS code and the INS2D code. The agreement with experimental 
data is fairly good for t/T < 2.0 for the two cases, but declines in quality for t/T > 2.0. 

Rosenfeld suggested that this disparity occurs as a result of the interpretation of the 
experimental results (ref. 14). The computed streamline traces represent instantaneous 
streamlines. The experimental streamline cannot be visualized instantaneously due to the 
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exposure time of the camera taking the streamline data. The estimation of actual physical 
parameters may affect the comparison of experimental and computed results. The effect 
of the error tolerances on the experimental results is unknown. 
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Another area of comparison is the computing time needed to obtain a periodic so- 
lution. Rosenfeld (ref. 16) reports the INS3D-FS code requires from 0.3 — 1.0 x 10 -3 
CRAY-XMP CPU sec/mesh-point/time-step, depending on the number of iterations re- 
quired for the convergence of the Poisson solver. To obtain a periodic solution, the INS2D 
code required 0.24 x 10 -3 CRAY-YMP CPU sec/mesh-point/time-step for Case I, and 
0.32 x 10~ 3 CRAY-YMP CPU sec/mesh-point/time-step for Case II. These cases required 
about 12 sub-iterations per physical time step. For the purpose of comparison, these cases 
were computed using the same size physical time steps and the exact same grids that were 
used in the study by Rosenfeld (ref. 14). Overall, the two different algorithms require the 
same amount of computing time to obtain a periodic solution. 


CONCLUSIONS 


A general study to compare and contrast two different algorithms for unsteady, in- 
compressible Navier-Stokes equations was performed. This was accomplished by using 
the INS2D code to generate computational results for comparison with the computational 
results of Rosenfeld (ref. 14) and the experimental results of Park (ref. 15). The compu- 
tational work done by the authors included several parameter studies and the generation 
of converged, periodic solutions. First, a grid resolution study was performed and an 
appropriate grid for the constricted channel configuration which sufficiently captured the 
important flow physics was generated. It was necessary to choose a finely spaced grid with 
289 x 61 grid points in order to compute significant flow structure details which were not 
adequately resolved by the coarser grids. In addition, the effect of the artificial compress- 
ibility coefficient, 0, on the speed of convergence and the smoothness of the solution was 
examined. For this configuration, a large value of 0 increased the speed of convergence but 
added “spikes” to the solution; as the value of 0 was increased, the magnitude of the spikes 
in the solutions increased. Several ideas were generated to reduce or eliminate the spikes 
in the solution. Grid refinement was performed by fitting a cubic spline to the streamwise 
grid lines in the area near the constriction. This grid refinement significantly reduced the 
magnitude of the spikes in the solution. Next, a study was performed to determine the 
number of periods the solution must be run out in order to achieve periodicity. After four 
pulsatile inflow periods, the flow physics in the solution became periodic. The effect of 
varying the maximum allowable tolerance for the divergence of velocity on solution accu- 
racy was examined. Several cases were run with different values of this parameter, £ con <, 
and the effect of these different values on accuracy and the amount of CPU time needed 
to generate a solution was determined. As e CO nt was reduced, the solution became more 
accurate, but also more computationally expensive to compute. In this study, a crossover 
point occurred at which further reducing the value of £ CO nt did not increase the accuracy 
of the solution, but was computationally more expensive. Finally, two converged, periodic 
solutions were generated and data was saved at each time step in the period for comparison 
with the computational results of Rosenfeld (ref. 14) and the experimental results of Park 
(ref. 15). 
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In comparing the computational results generated by Rosenfeld (ref. 14), using the 
INS3D-FS code, and the current work, using the INS2D code, it appears that the solutions 
generated by the INS3D-FS code contain stronger vorticity than those generated by the 
INS2D code. One area where this is evident is in the skin friction data. The skin riction 
data shown in figures 12(a)-(d) shows similar locations of the fluctuations m skin friction 
throughout the period. However, the magnitude of the fluctuations is generally larger in 
the solutions generated by Rosenfeld (ref. 14) using the INS3D-FS code, which correlates to 
stronger vorticity in these solutions. Due to a lack of detailed 

it is difficult to determine whether the solutions generated by INS3D-FS or INS2D are 
more correct in terms of flow physics. This dissimilarity is most likely a result t of the 
differencing schemes used by the two methods for the convective terms. The INS3D-tb 
code uses second-order central differencing on a staggered grid which does not add any 
artificial dissipation. The convective terms in the INS2D code are differenced using a 
third-order upwind scheme which does add dissipation. This additional dissipation may 
account for the difference in vorticity found in the two solutions in the region behind the 

constriction. 


The differencing schemes used for the convective terms comprise only one exam- 
ple of where the two algorithms are distinct, as the overall approaches used in the two 
methods are vastly different. The common thread in the two algorithms is the starting 
point of using the incompressible Navier Stokes equations. Starting from these equations 
the two methods follow radically different paths. The INS3D-FS code uses an Integra 
formulation of the equations, while the INS2D code uses a differential formulation. An- 
other difference in the two methods is the degree of coupling of these equations. In the 
current method, the INS2D code, the equations are fully coupled, whereas the equations 
used in the INS3D-FS code are only partially coupled through the pressure equation. The 
divergence of velocity criteria is satisfied in two completely different ways by the two meth- 
ods The use of a staggered grid approach by the INS3D-FS code, versus a non-staggere 
grid approach by the INS2D code, included with the use of different implicit schemes and 
dissimilar differencing schemes, all contribute to the distinctness of the two methods. How- 
ever, despite all the dissimilarities found in the two methods, the two algorithms generate 
very similar results, as shown in figures 15 and 16, which illustrate the location of the 
B-vortex and compare the computed results with experimental data. The two algorithms 
also generate very similiar streamline plots throughout the period. Overall, the agreement 
between the computed results of Rosenfeld (ref. 14) and the current work is excellent, 
and the two algorithms require similar amounts of computing time to obtain a periodic 

solution. 


FUTURE WORK 


One area which could be examined is the degree of robustness of the two algorithms. 
Future work could include checking the robustness of the two codes for larger physical time 
steps. This could be accomplished by continually increasing the size of the time steps by 
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a factor of two until the codes suffer a loss in stability. Another area which merits further 
investigation is grid sensitivity. As the two algorithms employ greatly contrasting spatial 
differencing schemes, examining the effects of skewness, discontinuities, and singularities 
in different grid topologies could provide interesting insights into algorithm characteristics. 
An additional way of examining the robustness of the two codes would be to continually in- 
crease the Reynolds number while monitoring the physics as well as the numerical stability 
of the algorithms. 
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